Effects of cold dark matter decoupling 
and pair annihilation on cosmological perturbations 
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Weakly interacting massive particles are part of the lepton-photon plasma in the early universe 
until kinetic decoupling, after which time the particles behave like a collisionless gas with nonzero 
temperature. The Boltzmann equation for WIMP-lepton collisions is reduced to a Fokker-Planck 
equation for the evolution of the WIMP distribution including scalar density perturbations. This 
equation and the Einstein and fluid equations for the plasma are solved numerically including the 
acoustic oscillations of the plasma before and during kinetic decoupling, the frictional damping occur- 
ring during kinetic decoupling, and the free-streaming damping occurring afterwards and throughout 
the radiation-dominated era. An excellent approximation reduces the solution to quadratures for the 
cold dark matter density and velocity perturbations. The subsequent evolution is followed through 
electron pair annihilation and the radiation-matter transition; analytic solutions are provided for 
both large and small scales. For a 100 GeV WIMP with bino-type interactions, kinetic decoupling 
occurs at a temperature Td = 23 MeV. The transfer function in the matter-dominated era leads to 
an abundance of small cold dark matter halos; with a smooth window function the Press-Schechter 
mass distribution is dn/dhiM oc M" 1/3 for M < 10~ 4 (T d /10 MeV)" 3 M . 

PACS numbers: 95.35.+d, 95.30.Cq, 98.80.Cq 
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I. INTRODUCTION 

Weakly interacting massive particles (WIMPs) are 
perhaps the leading candidate for the cold dark mat- 
ter (CDM) making up most of the nonrelativistic mass 
density of the universe today pj. Although candidate 
WIMPs are 10 to 1000 times more massive than nucle- 
ons and have no electromagnetic or color charges, their 
cosmic histories share many parallels with nucleons. Af- 
ter their abundances froze out at ~ 10 -10 the number 
density of photons, both WIMPs and nucleons remained 
thermally coupled to the plasma by elastic scattering 
with abundant relativistic particles. Acoustic oscillations 
in the relativistic plasma imprinted oscillations on both 
WIMPs and nucleons. Eventually the plasma released its 
grip on both types of particles. For WIMPs, this event 
is called kinetic decoupling; for nucleons, recombination. 

During kinetic decoupling, friction between the WIMP 
gas and relativistic plasma led to Silk damping of small- 
scale waves similar to what happened much later for 
atomic matter at recombination. After the respective de- 
coupling periods ended, pressure forces (and in the case 
of WIMPs, shear stress) inhibited gravitational instabil- 
ity on small scales. Still later, both WIMPs and nucleons 
played a major role in galaxy formation. 

There are, of course, significant differences between the 
cosmic evolution of WIMPs and nucleons. Most evident 
are the quantitative differences: because their interac- 
tions are so weak, WIMPs decoupled from the plasma 
less than one second after the big bang. Consequently the 
WIMP acoustic oscillations appear only on a length scale 
vastly smaller (~ parsec scales) than the baryon acoustic 
oscillations. The physics of nucleon decoupling, as im- 
printed in the galaxy distribution^] and in the cosmic 
microwave background radiation [3j provides a powerful 



probe of cosmic parameters and inflationary cosmology. 
If it were possible to similarly measure fluctuations on 
the scale of WIMP acoustic oscillations, we would have 
a dramatic consistency test of the cosmological model as 
well as an astrophysical measurement of WIMP proper- 
ties. 

One way to constrain the parameters of cold dark mat- 
ter decoupling is to measure the mass function of the 
smallest dark matter clumps today 13- E3| • Such clumps 
would be far too diffuse to host observable concentra- 
tions of atomic matter. However, they might be observ- 
able through the products of the very rare WIMP- WIMP 
annihilations taking place in the cores of these objects. 
Diemand et al. 7] proposed that numerous Earth-mass 
clumps might survive to the present day and provide a 
detectable gamma-ray signal. The mass and abundance 
of these clumps depends on cosmic fluctuation evolution 
during and after kinetic decoupling. 

WIMPS and nucleons also differ in a qualitative man- 
ner which has important consequences for the evolution 
of fluctuations through kinetic decoupling. After recom- 
bination, elastic scattering is rapid enough for atoms (and 
the residual free electrons) to behave as a nearly perfect 
gas on cosmological scales. Baryons behave like a fluid. 
WIMPs, however, collide too infrequently to thermalize 
after kinetic decoupling. WIMPs behave like a collision- 
less gas. Different approximations to the evolution of this 
collisionless gas have led to different results for the small- 
scale transfer function of CDM fluctuations d, H, 13 ■ 

In the present paper, the transfer functions for CDM 
fluctuations are calculated starting from the full Boltz- 
mann equation describing elastic scattering between 
WIMPs and the relativistic leptons present before neu- 
trino decoupling and electron-positron pair annihilation. 
Because the momentum transfer per collision with non- 
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relativistic WIMPs is small, the Boltzmann equation tion, 
reduces to the Fokker- Planck equation describing dif- 
fusion in velocity space caused by elastic scattering, 
combined with advection and gravitational forces. The 
Fokker-Planck equation fully describes kinetic decoupling 
and the evolution of perturbations of any length scale 
without approximating the WIMPs either as a perfect 
fluid or fully collisionless gas. Although the solution 
of the perturbed Fokker-Planck equation is more diffi- 
cult than the solution of coupled fluid and collisionless 
Boltzmann equations, it is both numerically and analyt- 
ically tractable (with an excellent approximation) in the 
present case. 

After kinetic decoupling, two additional events have an 
effect on the CDM transfer function. The first is electron- 
positron pair annihilation, which changes the equation 
of state of the plasma thereby modifying the evolution 
of fluctuations. Although the effects are small, they can 
be analytically calculated. The more important event 
is the transition from a radiation- to matter-dominated 
universe occurring about 10 5 years after the big bang. If 
the photons and neutrinos are treated as fluids, it is pos- 
sible to get analytic results for the linear evolution all the 
way to low redshift which are accurate to a few percent. where 
With these results in place, using standard techniques it 
is straightforward to estimate the mass function of CDM 
clumps at high redshift. 
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where M. is the Lorentz-invariant scattering amplitude, 
fi x = / x (pi) and similarly for the other distribution 
functions, Ei is the energy of particle i, and /l = 
(2irh) 3 /l/2 is the occupation number. Pauli blocking 
must be included for the leptons but not for the neu- 
tralinos since the latter have a low density after chemical 
decoupling. Equation Q is relativistically covariant but 
gives only the effects of collisions; the effects of grav- 
itational perturbations will be added later. Assuming 
effectively massless leptons, the matrix element for slep- 
ton exchange is given in Appendix A of and may be 
written 
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II. EVOLUTION OF WIMP PERTURBATIONS 
THROUGH KINETIC DECOUPLING 



Weakly interacting dark matter particles are described 
by their phase space density, which obeys the Boltz- 
mann equation governing transport by collisions with lep- 
tons (during kinetic decoupling these are just electrons, 
positrons, and neutrinos). For definiteness we will take 
the WIMP to be the lightest neutralino x°j however the 
results are easily applied to other WIMP candidates by 
modifying the scattering matrix element below. 

Let / x (p) and /l(p) be the proper phase space den- 
sities of neutralinos and ultrarelativistic leptons, respec- 
tively, where p is the proper three-momentum in a local 
orthonormal frame. (The spacetime coordinates x and t 
are suppressed for brevity.) The phase space densities are 
normalized so that J /(p) d 3 p is the spatial number den- 
sity, summed over spin states (we assumed unpolarized 
spins). One distribution function suffices for the rela- 
tivistic leptons because electroweak interactions maintain 
local thermal equilibrium at a temperature Tl. Elastic 
scattering of neutralinos and leptons cause the neutralino 
distribution to evolve according to the Boltzmann equa- 



ls a dimensionless constant depending on the relevant 
particle masses and couplings (Gf is the Fermi constant, 
bh and cl are left and right chiral vertices, and my/, 
77i£, and m x are respectively the masses of the W boson, 
the slepton, and the neutralino; Gprn^ = 0.0754). Here 
we assume following Ref. £| that the neutralino is a pure 
bino. Additionally, pcm is the momentum in the center of 
momentum frame, and (using metric signature — h H — h) 
tu = —{P2 — Pi) 2 — — 2p CM {\ — cos^cm) is one of the 
Mandelstam variables. 

For p CM ~ T L < m x , p C M = m x e[l - e + §e 2 + 0(e 3 )] 
where e = —p% ■p^jm^. In the lab frame, working to first 
order in Ti/m x , assuming p\ ~ x fm x Ti J ^ the collision 
kinematics gives 

1 - cos 6>cm =1 (n 2 + n 4 ) ■ (p 1 +p 2 ) 
1 — ri2 • 114 m x 

+ ( — ) (M?2 + + Ml4 - 1) . (4) 

where /zy = n.j • and rij = p;/p;. Another useful 
relation follows from energy conservation, 
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valid again to first order in T^/m x . 

Frequent collisions among the leptons maintain ther- 
mal equilibrium. Assuming negligible chemical potential, 
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for each species of massless lepton we have 



\fL(p)]- 1 = l + exp 



(l-n-v L ) 



(6) 



where is the (very small) local lepton fluid velocity 
due to cosmological perturbations. It is easy to check 
that an equilibrium solution of is then the Maxwell- 
Boltzmann distribution with mean velocity vx. 

For Tl <C m x the fe x term may be Taylor-expanded 
in (JJJ. After a lengthy calculation using l@J|-|JEJ, one 
obtains p^d^f = m x (df /dt) c (dropping the subscript on 
f x ) where the Boltzmann collision integral becomes the 
Fokker-Planck operator, 
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dp 



(P ~ m x w L )f + m x T L 



dp 



where 
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is a rate coefficient (in units where % = c = 1). Our exact 
result for the rate coefficient is larger by a factor 9.9 than 
the estimate obtained from Eqs. (9) and (12) of Ref. Q| 
and by a factor 3.4 than Eq. (17) of Ref. 5}. The rate is 
greater than the simple estimates made in previous work 
because of the details of the kinematics and the near- 
cancelation of forward and inverse rates in A larger 
rate coefficient leads to a lower temperature for kinetic 
decoupling than previous estimates. 

If we neglect spatial inhomogeneities, the unperturbed 
phase space density /o(?, r) depends on both comoving 
momentum q — ap and conformal time r according to 
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f . 2 rr ®fo 
q/o + a m x T L-g^ 



(9) 



Amazingly, for any time-dependence of 7, a, and Tl 
an exact solution to this Fokker-Planck equation is the 
Maxwell-Boltzmann distribution 



fo{q,r) 



exp(-g 2 /2cr 2 ) 
(2^)3/2 



o-g(r) = a(r), m x T x (T) , 



(10) 

where the WIMP temperature T x follows from integrat- 
ing 



dln(a 2 T x 
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During adiabatic evolution in the early universe, Tl oc 
a^ 1 oc t^ 1 and the WIMP proper temperature is then 
given in terms of the incomplete Gamma function by 



T x (r)=T iS VV r( 3 s) 



1 
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V aT = 2H 



(12) 



Equation Hl()|) may be multiplied by any constant, al- 
lowing the comoving number density of WIMPs to be 
normalized to its value after freeze-out. 



Familiarity with Brownian motion makes it seem nat- 
ural that the solution to the Fokker-Planck equation is 
a Maxwell-Boltzmann distribution. However, the lep- 
ton temperature Tl and the momentum transfer rate 7 
are falling with time and WIMP-WIMP elastic scatter- 
ing is far too slow to thermalize the distribution. Even 
so, collisions with the leptons maintain the WIMPs in 
a thermal distribution with a temperature that deviates 
increasingly from the lepton temperature throughout ki- 
netic decoupling. Once kinetic decoupling is complete 
the WIMP momenta redshift as p oc a -1 preserving the 
Maxwell-Boltzmann distribution with T x oc a~ 2 . 

Long before kinetic decoupling (s ^> 1), T x /Tl = 
1 - l/(4s) + 0(s~ 2 ). After kinetic decoupling (s < 1), 

x/T L ^ r(§) 

pling time by s = 1 yields 



T X /T L -> r(|)s 1 / 4 oc a- 1 . Defining the kinetic decern- 



T d ^ Tl(s = 1) = 1.430 C-^g^i^L 



1/4 



m P i 



7.650 C~ 



For typical supersymmetry masses, kinetic decoupling 
occurs after muon annihilation when the only abundant 
leptons are electron pairs and neutrinos, for which (with 
sin 2 9 W = 0.223) + c£) = f tan 4 9 W = 0.268 and 

9cS = 43/4. With m x = 100 GeV and m L = 200 GeV, 
C = 0.0433 yielding T d = 22.6 MeV. Profumo et al. [13 
show that Td can span the range from a few MeV to a 
few GeV for reasonable WIMP models. Here we take the 
supersymmctric bino as a candidate but will show how 
numerical results scale with T^. 

Next we examine the effect of density, velocity, and 
gravitational potential fluctuations during kinetic decou- 
pling. The perturbed phase space density is /(x, q, r), 
where x are comoving spatial coordinates, q = ap are 
the conjugate momenta, and r is conformal time. The 
perturbed line element in conformal Newtonian gauge is 
written ds 2 = a 2 [-(l + 2$)dr 2 + (1 - 2*)dx • dx]. In- 
cluding the effects of the metric perturbations $ and ^, 
the Boltzmann equation becomes 

9f df (■ d<§>\ df fdf\ 



where v = q/e is the proper velocity measured by a 
comoving observer, e = (q 2 + a 2 ?™ 2 ) 1 / 2 is the comov- 
ing energy, and q = gn. With comoving variables, the 
Fokker-Planck operator becomes 
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df 

(q-qL)/ + a 2 m x T L ^ 



(15) 



where q^ = am x vi = —am x \JuL- 

Now we linearize 111 for first-order perturbations of 
the lepton fluid by writing T L -> T l {t) + 7\(x, r). The 
fields Ti(x, r), $(x, t), ^(x, t), and the lepton velocity 
potential ul(x, r) are treated as first-order quantities. 
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Assuming v 2 <C 1 and performing a spatial Fourier trans- 
form, we obtain the linear perturbation equation 



1 

To 



df 
dr 



«(k- v)/ - 7 a ^Fp/ 



■ Q «k ■ q , . ( q 2 

= — $ + jau L ) + 7 a — 

a 2 aT x \a- 



where 
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is given by k = k/fc . Thus we may write 
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(27ra 2 m x T L ) 3 / 2 
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x J2 (-i)\2l + l)Sni(y)Pi(k-n)f nl (k,T) , 



,1=0 



S nl EE ^/2 L ( i+ l/2) (2/) . 



(22) 



Before writing the perturbation equations, let us ex- 
amine the unperturbed solution, / = fo(q, r)5 3 (k). The 
exponential factors differ in I jlOjl and l)22[l. implying that 
the Laguerre expansion includes more than one term. In- 
deed, one finds 



After WIMP freeze-out, the leptons dominate the gravi- 
tational potentials so that ^ , ul, and T\ are functions 
of (k, r) given by the solution for a perfect relativistic 
fluid. Equation (|16[) generalizes the collisionless Boltz- 
mann equation of Ref. to include the effects of dark 
matter collisions with relativistic leptons. 

To integrate l(TB|l for the phase space density / we will 
expand the momentum dependence in eigenfunctions of 
the Fokker-Planck operator Lfp: 

L-evfynim = -(2n + l)4> n i m , (18a) 
Ki m = e-y /2 4 i+1/2) (2/)>Wn) • (18b) 

Here y = q 2 / (2a 2 m x TL) and L n a \y) is a generalized La- 
guerre polynomial, also known as a Sonine polynomial. 
It is defined by the following series expansion: 



y^ r(n + Q+l) (-X) k 

^T{k + a + l) k\(n-k)\ 



(19) 



Generalized Laguerre polynomials have orthonormality 
relation 



Jo T(n + a + l) 
and completeness relation 



L^\x)L^\x)dx = 5 nnl (20) 



E 



n\ (xx') 



Ao/2 p -(x+x')/2 



T(n- 



1) 



L£\x)LM( X ') = 6(x-x') 



(21) 

We will expand the phase space density /(k, q, r) in 
the complete set 4> n lm- However, it is unnecessary to 
include all (nlm) in this expansion. Prior to kinetic de- 
coupling the Fokker-Planck operator Lpp rapidly damps 
all terms except </>ooo and those terms that are induced 
by the right-hand side of i|16|) . The rotational symmetry 
of this equation implies that only the m = terms are 
induced, where the polar axis for the spherical harmonics 



/ = /o(<?,t) 



fnl — <>U 



1 ^ 



(23) 



Prior to kinetic decoupling, when collisions maintain 
T x = Tl, /oo = 1 and the other coefficients vanish. Af- 
ter kinetic decoupling, /„o — * 1 for all n. The Laguerre 
expansion must be carried to high order in order to con- 
vert e~ v to exp(— q 2 /2a 2 ). Similarly, we should expect 
the perturbed phase space density also to require many 
terms in n after kinetic decoupling. 

Substituting l|22|l into fTSjl and using orthonormality 
and several recurrence relations for the generalized La- 
guerre polynomials, we obtain a system of coupled ordi- 
nary differential equations for the evolution of the per- 
turbed phase space density, 

fni + (2n + 0(7« + R)f n i - 2nRf n -i,i 

+k ^{ (mt)[( n + l + D/«,i+l " »/n-l,I+l] 

+ 2fp[(fn+l,l-l ~~ /n,i-l)| 

= 5 l0 [WA n - 2(* + 7 aTi/T x )fl n ] 
+hi - \ I -7^- + iau L )A n , 



where 



R 

Au{t) 
B n (r) 



To 



(24) 



(25) 



The density perturbation is <5(k, r) = /oo for k ^ 0. 

The source terms for Eqs. (|24f) are provided by the rel- 
ativistic plasma. Their time-dependence changes during 



1 Collisional damping before neutrino decoupling similarly justifies 
the neglect of the in ^ modes for neutrinos in Ref. Illl . 
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lepton pair annihilation and neutrino decoupling. For 
reasonable parameters, neutralino kinetic decoupling oc- 
curs after muon annihilation but before electron annihi- 
lation and neutrino decoupling. During this era, the isen- 
tropic mode of perturbations has time-dependence given 
by the relativistic perfect-fluid transfer functions 



$ = * 



Hul = - 



3 . 
-^3 (sin ( 

3 

29 



1 - 



9 cos ( 
2 



3 



-<P-Hu L , 



(26a) 
(26b) 

(26c) 



where TL = a/a, 9 = kr/y/3, and the transfer functions 
are normalized so that <f> = — 1 for = 2 . The actual 
perturbations are obtained by multiplying the transfer 
functions with the scalar field (k) that gives the spatial 
dependence of the initial (inflationary) curvature fluctu- 
ations. As the inflationary curvature perturbations are 
well known (a Gaussian random field with nearly scale- 
invariant spectral density P cx fc -3 ), here we work with 
transfer functions. 

Initial conditions for Eqs. I|24|l are obtained by exam- 
ining the solutions for s = \^o,t = (t/t^) -4 ^> 1. Isen- 
tropic initial conditions have /oo = S = 5^. The only 
significantly perturbed components of / for the strongly 
coupled plasma are the thermal equilibrium values 



foi 



2™, 



TV. 



UL , /lO — = -n^L ■ 

(27) 

All other components /„; are kept small by rapid WIMP- 
lcpton collisions. Equations (|24|l with n < 1400 and 
I < 15 were integrated to r = 72rd using the standard ex- 
plicit ordinary differential equation solver DVERK. Con- 
vergence testing showed that higher-order terms in the 
Laguerre expansion were negligible. 

Figurenshows the results for the density transfer func- 
tion expressed using the gauge-invariant variable v, de- 
fined as the CDM number density perturbation in a syn- 
chronous gauge for which the mean CDM velocity van- 
ishes in the coordinate frame (for a nonrelativistic parti- 
cle, v equals Bardeen's variable e m ). This variable, which 
is used by CMBFAST llj, is related to the conformal 
Newtonian gauge variables by 



v — 8 + ZTLu 



(28) 



For wavelengths longer than the radiation acoustic length 
r/V3, v = 5 L + 3TCu L ^(kr) 2 . For wavelengths shorter 
than the radiation acoustic length at r but longer than 
the acoustic length at (i.e., 10t _1 < k < r7 ), the 



2 Here, 8l = &Pl/{pl + Ph) is the number density perturbation 
in conformal Newtonian gauge; the energy density perturbation 




for the relativistic leptons is Sp^ /pL 



4 x 



FIG. 1: CDM density transfer function versus wavenumber 
at conformal time r = 72-Tcj. The three oscillating curves as- 
sume that kinetic decoupling occurred at Ta for three different 
values of the radiation temperature Td relative to the CDM 
particle mass m x ; the amplitude of the oscillations decreases 
with increasing Td/m x . The upper, monotonic curve assumes 
that the CDM is always collisionless and was never coupled to 
the radiation. The other non-oscillating curve shows a crude 
model of kinetic decoupling described by a Gaussian cutoff. 



acoustic oscillations of the gravitational potential aver- 
age out leading to a suppression of growth induced in the 
CDM. For these intermediate wavelengths the transfer 
function is a logarithmic function of wavenumber. If the 
dark matter were completely non-interacting, this loga- 
rithm would continue to arbitrarily high wavenumbers, 
as illustrated by the monotonically rising curve in Fig. ^ 

Because WIMP dark matter was collisionally coupled 
to the relativistic lepton plasma at early times, the CDM 
transfer function in Fig. ^ shows remnant damped acous- 
tic oscillations at short wavelengths. For comparison 
a Gaussian transfer function is shown, with no oscilla- 
tions. In this model, the effects of kinetic decoupling 
are described by multiplying the transfer function for 
the completely non-interacting case by exp(— k 2 /2a 2 ). 
As we will see, a simple model of free streaming pre- 
diets that KTd)- 1 « [r(|)T d /m x ] 1 /2[l n ( r / T(i )] dur- 
ing the radiation-dominated era. For r = 727^ and 
(2T d /m x ) 1 / 2 = 0.01 this model would predict <r k T d w 30 
whereas the curve shown in Fig. ^ has o^Td — 18. Free- 
streaming does not give a good approximation to the 
actual transfer function. 

The exact transfer functions shown in Fig. ^ decrease 
less rapidly with wavenumber than the approximate 
transfer functions of Ref. which were computed using 
a fluid approximation followed by free-streaming. The 
following section reviews the free-streaming solution and 
then develops a more accurate approximation based on 
moments of the exact Fokker-Planck equation. 
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III. APPROXIMATE DESCRIPTIONS OF 
PERTURBATION EVOLUTION THROUGH 
KINETIC DECOUPLING 

In this section we consider two different approxima- 
tions which provide analytical insight to the numerical 
solution of the Fokker-Planck equation. The CDM be- 
haves at early times like a fluid and at late times like 
a free-streaming collisionless gas, and in these limits we 
can develop useful analytical approximations. 



A. Free-streaming model 

For r ^> Td, the terms proportional to 7a may be 
dropped in The differential equation may then be 

integrated to give /(k, q, r) in terms of the initial value 
/(k, q, r*) for any 3> Td- Integrating over momenta 
gives the conformal Newtonian gauge density perturba- 
tion, 

5(k,r) = yd 3 g e- ,;ck q /(k,q,r,)+^ dr' e- M ' I 2 

x (3-Af)*(k,T') + £:Vln(^) $(k,r') , (29) 

where we have assumed evolution in the radiation- 
dominated era with 







a*m x 



In 



(-) , M'^^{kr'f^i-) ■ (30) 



For large spatial frequencies, fcr 3> 1, the gravitational 
potentials — which are dominated by relativistic parti- 
cles — oscillate rapidly leading to a small net integral 
contribution; ignoring this, <5(k, r) is a momentum-space 
Fourier transform of the distribution function at the ini- 
tial time r*. 

Obtaining the exact solution still requires numerical 
integration of (|16() through kinetic decoupling, or equiv- 
alently, integrating the system of equations H24JI . How- 
ever, we can get an idea of the effects of free-streaming 
by making an instantaneous decoupling approximation, 
treating the CDM as a fluid with a Maxwell-Boltzmann 
velocity distribution for r < t* and by l|29(l for t > r*, 
as was done in Ref. [9j. Then, whether or not the CDM 
is strongly coupled to the radiation, the perturbed dis- 
tribution function is fully specified by the perturbations 
of density, temperature, and velocity v = — iku, 



Sf = fo(q,r) 



^ + 1(^-3 



ST X iuk ■ q 



(31) 



P " \ u q / ^X "^X 

Carrying out the Fourier transform in l|29|) now yields 



S(r) 



= e- M ' 2 



<S-ttS^) - tfV.ufc) In ( — 



2 T v 



r* 



where wavenumber arguments have been dropped for 
brevity, and 



jU-2 _ 



rlT d 



(33) 

If t* ~ rd, then ST X /T X m |<5; if t* ^> and the fluid is 
approximated as being adiabatic, ST X /T X — |<5. 

Equation l|32|) corrects errors in the definition of kf 
of Ref. |9j and adds the term proportional to 5T X /T X to 
their Eq. (20). This new term arises from the tempera- 
ture perturbation of the CDM fluid, which modifies the 
distribution function and, through free-streaming, mod- 
ifies the density for t > t*. In particular, if the CDM 
perturbations are approximately adiabatic, the tempera- 
ture perturbation causes the transfer function to decrease 
more rapidly with k. 

The free-streaming solution predicts a Gaussian cutoff 
of the transfer function, S oc exp(— k 2 /2a 2 ), with cutoff 
distance equal to the free-streaming distance, 




dr = 



r(|) 



-,1/2 



dr 
a/a d 



(34) 

During the radiation-dominated era, when a oc r, the 
free-streaming length grows logarithmically, but it satu- 



rates in the matter-dominated era when 



At first 



(32) 



glance, Fig. ^ appears to qualitatively support the free- 
streaming model of a Gaussian cutoff. However, the free- 
streaming model predicts no damping for a super-heavy 
particle with Td/m x — > 0, while Fig. ^ shows that even 
in this case there is damping. This damping arises from 
friction between the lepton and CDM gases during de- 
coupling (Silk damping). This friction can be accounted 
for by treating the CDM as an imperfect fluid. 



B. Imperfect fluid model 

To better describe an extended period of decoupling 
while allowing for small deviations from a Maxwell- 
Boltzmann distribution, we consider the evolution of the 
lowest order moments of the distribution function. Wc 
work to first order in perturbed quantities and normalize 
the unperturbed distribution function to J fod 3 q = 1. 
The perturbations for k ^ then define the lowest-order 
moments 



fd 3 q = 5, 
fvj d 3 q = —ikjU , 
fviVj d 3 q = (c 2 x 8 + a)Sij - l(kikj - i/c 2 %)7r , (35) 



where the density perturbation S, velocity potential u, 
shear stress potential 7r, and entropy perturbation a are 
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related to our expansion coefficients by 



8 = /oo , ku= a ^— ^ fa , k 2 ir = — -f 02 , 



» 2m x 
c /oo J 10 



m. 



(36) 



The effective sound speed squared of the CDM fluid is 



T x ( x _ld\nT x 



3 din a 



(37) 



which differs from the thermal speed squared Tl/to x 
appearing in 1)24(1, This difference arises because the 
Laguerre expansion uses eigenfunctions of the Fokker- 
Planck operator which depends on the relativistic lepton 
temperature Tl rather than the WIMP temperature T x . 
After kinetic decoupling, T x /Tl drops and the higher- 
order expansion coefficients /„; will increase to compen- 
sate for this difference, as we already found happening 
with the unperturbed distribution function in l|23(l. 

The variables in l|3bj) describe fluctuations of an im- 
perfect fluid. The reader may wonder how a single 
component fluid can have an entropy perturbation. A 
weakly imperfect fluid is described by an equation of state 
p = p(p, S) where S is the entropy which may vary in 
space and time. However, the WIMP gas is more com- 
plicated because it becomes fully collisionless after ki- 
netic decoupling; it may be regarded as a superposition 
of many non-interacting ideal gases. 

The time evolution of the imperfect fluid variables fol- 
lows from Eqs. J2IJ 3 : 



ii - 



k 2 u = 
Hu = 

-2Ha 



3* , 
* + c 2 x 8 

5Tl 
3 m x 



- c~ | fc u 



(38a) 

k 2 n — ja(u — ul) , (38b) 

3/2 

kfu 



4\m x J 



-27a 






a — 






m x 


4 Tl_ 




/27V 




u — 




3 m x 






—2jan . 





:W I - — 

3 m v 



3/2 



~ c 'x I 5 



2I_io3 

4 k 



(38c) 



fn 
k 



(38d) 



These equations are similar to those of an imperfect gas 
coupled to the lepton fluid. However, they differ sig- 
nificantly from the Navier-Stokes equations assumed in 
Ref. In place of a bulk viscosity term (k 2 u and a 



3 These fluid equations, like the Fokker-Planck equation from 
which they were derived, are correct only to leading order in 
T L /m x . The corrections to I38al and I38bl are ^ — » ^ - Ti<r 
and S — > v = 5 + 3H.U, respectively. 



shear viscosity term ^i]k 2 u where £ and r\ are the bulk 
and shear viscosity coefficients (divided by the mass den- 
sity), l(38bp has an entropy term a and a shear stress 
term —k 2 ir where a and ir are not proportional to u. 
The usual Chapman-Enskog expansion does not apply to 
our Fokker-Planck equation when the collision mean-free 
path becomes large. Moreover, because of the /n and /03 
terms, Eqs. (I38|l do not form a closed system. Nonethe- 
less, these equations are useful for providing insight and 
they will guide us to a very good approximation to the 
full numerical solution of the Fokker-Planck equation. 

Prior to kinetic decoupling, when the damping terms 
proportional to 7a are large, the solutions to (|3~5|l have 
5 = Sl, u = ul, a = n = 0, in agreement with l|27|l . 
Entropy and shear stress perturbations develop during 
decoupling as the CDM gas becomes collisionless. These 
in turn modify and damp the acoustic oscillations of the 
gas. 

It is instructive to solve the imperfect fluid equations 
with several different approximations, in order to deduce 
which physical effects are responsible for the features of 
the transfer functions shown in Fig.^ The most extreme 
approximation is to completely neglect the CDM temper- 
ature and coupling to other particles, so as to describe a 
perfect cold, collisionless gas. This approximation con- 
sists of setting c 2 , = c = 7r = 7 = retaining only the 
gravitational interaction between the CDM and the rel- 
ativistic plasma. In this case the exact solution prior to 
neutrino decoupling for isentropic initial conditions is 
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COS ( 



1 sin ( 



Ci(0) +ln0 



Hu = — r(sin0 — 



(39) 



where 9 = kr/VS, C = 0.5772 •• • is the Euler- 
Mascharoni constant, and Ci(9) is the cosine integral. 
This solution is shown by the monotonically increasing 
curve in Fig.^ We see that once the wavelength becomes 
smaller than the relativistic acoustic horizon and the po- 
tentials oscillate faster than the CDM can respond, the 
CDM density perturbation growth slows to logarithmic 
in time. This suppression is responsible for the turnover 
of the low-redshift CDM power spectrum from P(k) oc k 
at long wavelengths to P(k) oc fc -3 ln 2 (fc) at short wave- 
lengths. However, the approximation of a cold, collision- 
less fluid includes none of the physical effects of kinetic 
decoupling. 

The next simplest approximation is to treat the CDM 
as being cold (c x = a = ir = 0) but to include the fric- 
tion term — ^a{u — ul) in the velocity equation. This 
approximation is exact in the limit Td/m x — > hence it 
reproduces the case T ( j/m x = in Fig. ^ We can find 
the solution in the limit t > by first noting that the 
general solution of the second-order system 5+k 2 u = 3^f, 
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0.0 


1.0000 


1.0 


-0.1637 


2.0 


0.0252 


3.0 


0.0099 


4.0 


-0.0101 


0.1 


0.5744 


1.1 


-0.1562 


2.1 


0.0345 


3.1 


0.0039 


4.1 


-0.0082 


0.2 


0.3791 


1.2 


-0.1409 


2.2 


0.0403 


3.2 


-0.0013 


4.2 


-0.0062 


0.3 


0.2275 


1.3 


-0.1202 


2.3 


0.0427 


3.3 


-0.0056 


4.3 


-0.0042 


0.4 


0.1063 


1.4 


-0.0965 


2.4 


0.0424 


3.4 


-0.0090 


4.4 


-0.0022 


0.5 


0.0109 


1.5 


-0.0716 


2.5 


0.0397 


3.5 


-0.0113 


4.5 


-0.0004 


0.6 


-0.0610 


1.6 


-0.0471 


2.6 


0.0352 


3.6 


-0.0126 


4.6 


0.0011 


0.7 


-0.1196 


1.7 


-0.0244 


2.7 0.0295 


3.7 


-0.0130 


4.7 


0.0023 


0.8 


-0.1443 


1.8 


-0.0045 


2.8 


0.0231 


3.8 


-0.0126 


4.8 


0.0033 


0.9 


-0.1606 


1.9 


0.0122 


2.9 


0.0164 


3.9 


-0.0116 


4.9 


0.0039 



TABLE I: Residuals from the 
to fi(x), defined in Eqs. tUl 



steepest descent approximation 
and lEESl . 



u + TLu = $ is given by 
cos 8 + 8 sin 8 



i' 
9 

Uu 



ff 2 

{smO-fxt 



Ci{6) + h[\n6- ¥ 



8-- 



(40) 



where f± and ji are independent of r but may depend 
on k. Including the damping terms in i|38bjl promotes 
fi and to functions /i(fc,r) and fiik,T) obeying the 
differential equations 



fi + jafi 
h+ftlnO = 



" I cos 8 + —8sin( 



(41a) 
(41b) 



In the strongly-coupled limit t <C Td the solution must 
match Eqs. I|26(l . which gives 



h(k,r) 



cos 8 H — 8 sin 8 , 



Ci(8) 



cos8 + -8sin8 ) In 6 



■ cos 8 



(42) 



The exact solution of l|41|) satisfying these initial condi- 
tions is given by a pair of quadratures, 



fi(k,r) 
f 2 (k,T) 

where s : 



e( s - s '>/ 2 | cos 



0' + -0'sin0' hVrfr' 



dr' 

[ fl (k,r')-l] — + (l-f 1 )hx6 

T 



(43) 



\^aT = (T/Td)^ 4 , 8 = kr/\/3, and primed 
quantities are evaluated at r'. 

The late-time solution is dominated by /i, which for 
t Td becomes the following function of wavenumbcr, 

h{x) =Ke e- z{r ^ '-dz, x = I —4 1 



V2V3/ 



(44) 



N 

E 




FIG. 2: Contour used to evaluate the integral in Eq. 144^ . The 
quarter-circle is actually taken to have a much larger radius 
than shown so that its contributions to the contour integral 
vanish. The desired path for the integral is along the lower 
left curve that joins the quarter-circle. 



The complex function z{r) — xe~ ' (^r + 2r) where 
r itself is complex but cannot be used as the integration 
variable because of the essential singularity at r = 0. 
The contour used to evaluate /1 is shown in Fig. [3 Us- 
ing the steepest descent approximation to evaluate the 
contributions along the branch cut gives 



/ 47ra N 




5 , 


X 


cos 


5.?; 






sin 


T 





1/2 



exp 



5 /2tt 
--xcos i — 



3tt 

X COS [ <y9 — 



ip = — sin — 



(45) 



where f rcs (x) is a correction to the steepest descent ap- 
proximation, evaluated numerically and tabulated for in- 
terpolation in Table [I] For x <C 1, 1 — /1 goes to zero 
faster than x 4 , with /i(0.2) = 0.9990. 

Numerical integration yields an approximate fit to 
f2{x) for x ^ 1 similar to l|45l) : 



f 2 (x) « -3.45x 



0.75 



5 

-—X cos 
2 



eN.p — — ./'cos( — j (x sin (p — cos (p) , 

(46) 

Figure 01 shows the results for f\(x) and f2(x). Although 
fi has larger amplitude, /1 dominates the contribution 
to the CDM transfer function for t ^ Td- 

This calculation shows that, in the limit of large WIMP 
mass m x 3> Td, the free-streaming suppression of the 
transfer function is not Gaussian but instead is exponen- 
tial in x oc fc 4 / 5 . 
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(2T d /m x ) 1/2 = 0.02 



30 40 
kx d 



50 60 



FIG. 3: The auxiliary functions fi{x) and fz(x) (the curve 
with the higher amplitude of oscillation) appearing in Eq. 
lEOt . where r > r d and x = (fcrd/273) 475 . 



Allowing a nonzero Td/m x introduces a Gaussian sup- 



pression, as follows. Including the c x S term in l(38 
adds the term — /c 2 c^r^/9 to the right-hand side of (|41a|l . 
With the leading late-time behavior e/9 « /iln#, this 
gives an approximate Gaussian suppression, 



v(k, t) « exp 



IP 

2k 2 



V (k,T) 



where ^o(fc,r) is the solution for Td/r 



2 x (r')dr' , 
(47) 

= given by 

gUJ) and |g3J|. While g7Jl is similar to (J22Jl and lO, the 

Gaussian cutoff distance differs from the free-streaming 
distance l|34JI by a factor 

Equations l|4(J|) and (|47|l assume that the CDM behaves 
like a perfect fluid 4 . By integrating (|38|l numerically 
with the simplification fn = /03 = 0, we can include 
effects of nonzero shear stress it and entropy perturba- 
tion a. The results are shown in Fig. 0] in comparison 
with the exact solution from Fig. ^ While the effects 
of nonzero temperature, shear stress, and entropy per- 
turbations qualitatively reproduce the suppression of the 
transfer function, none of the fluid approximations gives 
a good match to the exact solution of the Fokker-Planck 
equation. However, an excellent fit (with maximum error 
about one percent of the oscillation amplitude) is given 
by a modification of (|47|l . 



k' 1 



exp ' } t/ °( fc ' r ) ' 



6T d 



(It 



5m x J Tt a/ad 



(48) 



4 For example, Eq. 1471 describes the small-scale damping of fluc- 
tuations after recombination in a baryon-dominated universe. 



FIG. 4: CDM transfer function and several approximations 
plotted versus wavenumber at conformal time r = 72rd. In 
descending amplitude of the second peak, the curves are (1) 
fluid approximation with n — a — 0, (2) exact solution of 
the Fokker-Planck equation, (3) imperfect fluid approxima- 
tion with shear stress but no entropy perturbation; (4) imper- 
fect fluid approximation with both shear stress and entropy 
perturbations. The plus signs superimposed on the exact so- 
lution curve are the solution with Td/m x — 0, multiplied by 
a Gaussian damping factor as described in the text. 



This approximation, with = 1.05^, is shown by the 
plus signs in Fig.QJ The coefficient in front of the integral 
defining k^ 1 was found numerically; the value | = 1.200 
is numerically correct to 0.1% or better and clearly differs 
from the coefficient T(|) = 1.2254- • • appearing in the 
free-streaming prediction of l|34|) as well as the perfect 



fluid prediction |r(|) = 2.0423 • • • of gTJ. 

The modified form of thermal damping suggests that 
on small scales the WIMP gas might be described as 
a thermal gas with ratio of specific heats 1.2/P(|) = 
0.979 • • • instead of |. However, numerical tests showed 
that no perfect gas equation of state can reproduce the 
exact solution of the Fokker-Planck equation as well as 
(|48|l . We will therefore use l|48(l to calculate the effects 
of free-streaming damping even though it is based on a 
numerical instead of an analytic solution. Note that Ij48(l 
is valid through pair annihilation and radiation-matter 
equality because the free-streaming distance is propor- 
tional to j yjT x jm x dr and T x oc (a/a d )~ 2 at all times 
after WIMP decoupling. 

The physics of WIMP decoupling is similar but not 
identical to the much later decoupling of atoms at a tem- 
perature of 0.25 eV. In both cases the decoupled particles 
bear the imprint of acoustic oscillations while they were 
coupled to a relativistic gas. In both cases the acous- 
tic oscillations are damped by friction during decoupling 
(Silk damping), although the damping is exponential in 
fc 4 / 5 for WIMPs as opposed to k 2 for atoms. In both 
cases short-wavelength fluctuations are damped further 
by thermal motions after decoupling. However, the last 
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(2T d /mJ = 0, 0.01,0.02 




FIG. 5: Real space CDM Green's function (Fourier transform 
of the transfer function) at conformal time r = 10 7 Td, ap- 
proximately at the end of the radiation-dominated era. An 
initial planar perturbation sends an acoustic wave through 
the relativistic plasma. This wave travels through the CDM 
until kinetic decoupling ends; thereafter the wave is frozen 
in place but grows logarithmically in amplitude. The three 
curves show the effect of diffusion with increasing CDM tem- 
perature. 



stage differs for the two gases because the atomic gas 
remains collisional (hence damping takes place for wave- 
lengths shorter than the Jeans length) while the CDM gas 
is collisionless and the relevant scale is the free-streaming 
length. As a result the CDM gas develops shear stress 
and entropy perturbations in i|38bl) that are not present 
for a collisional gas. Nonetheless these perturbations lead 
to free-streaming damping that is qualitatively similar to 
the Jeans damping of a collisional gas. 

The physics of perturbation evolution can be simpler 
to interpret in position space than in Fourier space |l3j| . 
Figure [3] shows the one-dimensional Green's function de- 
fined by 



(x,t) 



1 

2^ 



v(k, t) dk 



(49) 



of neutrino free-streaming are small and are already in- 
cluded in numerical codes such as CMBFAST and COS- 
MICS but pair annihilation is not. As we show next, 
pair annihilation also results in a minor modification of 
the CDM transfer function on small scales. 



IV. EVOLUTION THROUGH PAIR 
ANNIHILATION 

During electron-positron annihilation the equation of 
state changes slightly, modifying the evolution of pertur- 
bations. Neglecting the small amount of neutrino heating 
that takes place during pair annihilation, after CDM ki- 
netic decoupling the energy density and pressure of the 
photons, 3 flavors of neutrinos, and electron pairs is 



Pi = 3p 



15 



Pv = %Pu 



7^ 4 o 
40a 4 



2T 4 
P± = -f 



2T 4 



(50) 



where £ = m e /T 7 and T v $ = aT v = constant. The pair 
density and pressure are given in terms of the Fermi- 
Dirac integrals 



Jo 



t 2 yjx 2 + £ 2 dx 
x A dx 



^x 2 +£ 2 (eV^+? +1) 
Energy conservation for the photon-pair plasma gives 
dln£ A 



(51) 



where 



.4 
B 



1 



1 



dim 



30R 



D 



(52) 



7T 

30i? 



15 

2^ 
15 dR 



4 (P-R) 



(53) 



r 4 2TT 4 dlny 
The photon-pair plasma has equation of state parameter 



In real space the Green's function is essentially a wave 
packet of sound that started at x = and propagated 
until kinetic decoupling after which time it froze in place. 
The oscillations in fc-space arise solely because the wave 
packet in real space has a rapid change in slope at its 
trailing edge. In the absence of WIMP-lcpton coupling, 
i^ 1 ) would have a delta function singularity at x = 
and a compensating underdense tail at x > due to 
the gravitational perturbation caused by the outgoing 
acoustic wave in the relativistic plasma |13|. 

Until now we have assumed that the photon-lepton 
plasma is a perfect relativistic gas with p = |p. This ap- 
proximation breaks down after neutrino decoupling and 
during electron/positron pair annihilation. The effects 



(10/7T 4 )(P-i?) 



1 + (30/7r 4 )i? + (21/8)(T V /T j y 



and sound speed squared 



3d 



A(A - B) 



S[A+(21/8)(T,/T 7 ) 4 



(54) 



(55) 



These are plotted in Fig. |SJ Pair annihilation makes a 
10% dip in the equation of state. 

The relation between expansion factor and conformal 
time follows from integrating the Friedmann equation 



47T 3 G(aT 7 ) 4 



45 



-9eS 



(56) 
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FIG. 6: Equation of state parameter w (solid curve) and ef- 
fective sound speed squared cj, (dotted curve) through the 
period of electron pair annihilation. 




FIG. 7: Amplitude C (solid curve) and phase shift a (dashed 
curves, normalized by 8 a = hr a /vZ) for the gravitational 
potential <[/ of Eq. I59L for r ^> r a , plotted versus d a (di- 
mensionless wavenumber). Three sets of curves are shown, 
corresponding to the times r = (10 3 ,10 4 ,10 5 )r a . Before pair 
annihilation, C — 1 and a = for all k. 



where 



30 ( p 2 + pu + p± \ 21/T„r 60 
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R . 

Neglecting the small neutrino shear stress arising from 
free-streaming after neutrino decoupling at T ~ 2 MeV, 
the conformal Newtonian gauge gravitational potentials 
$ = '3/ obey the evolution equation 

$ + 3(1 + <*,)H& + 3(4 - w)H 2 <$> + k 2 c 2 w $ = . (58) 

Long after pair annihilation is completed, w — c 2 w = i 



and the solution is 
3C 

IF 



$(fc,r) 



[sin(0 + a) -6 cos(8 + a)} 



(59) 



where C (not to be confused with the collisional coupling 
constant of Eq. |3J) and a are independent of r but may 
depend on fc. Before pair annihilation, C = 1 and a = 0. 
During pair annihilation, C and a change but after pair 
annihilation they become independent of time. 

We define a characteristic conformal time for pair an- 
nihilation, T a , by scaling from the decoupling time: 



0.511 MeV 



(60) 



Pair annihilation imprints features on $(fc,r) at 
wavenumbers k ~ t" 1 . After pair annihilation is com- 
pleted, we expect C and a in (|59|l to depend only on fcr a . 
Numerical integration of (|58|l yields the results shown in 
Fig. [7] Pair annihilation leads to a small change in the 
amplitude C of gravo-acoustic oscillations and a char- 
acteristic phase shift a of order 9 a = fer a /v3. Wave- 
lengths much longer than the acoustic horizon length 
are unmodified. Waves that en- 
ter the Hubble length when the effective sound speed 
is reduced by pair annihilation (Fig.|SJ) are amplified be- 
cause with decreasing sound speed, pressure forces are 
less able to prevent gravitational growth. The diminished 
sound speed also changes the distance traveled by acous- 
tic waves, leading to a phase shift for wavelengths smaller 
than the acoustic horizon distance. The propagation of 
the acoustic horizon is evident in the propagating steps 
in a /8 a in Fig. [7| Although these steps would appear to 
prevent <I> from relaxing to (|59|) with time-independent C 
and a, for r ^ r a , 6 \6 a \ and one may set a = with 
an error in $ of 0(r a /r). Wavelengths much longer than 
the Hubble length during pair annihilation are essentially 
unmodified. 

Still shorter waves are described by the WKB solution 

of 



$(fc, t) cx 



p + p 



1/2 



COS 



cc' dr' for fcr > 1 



(61) 

The dip in the sound speed during pair annihilation leads 
to a phase shift a = k f^(c' w - 1/V3) dr' = -0.4056» a for 
9 a ^> 1 and r ^> r a . Pair annihilation effectively resets 
the starting time for acoustic oscillations from t = to 
t = 0.405r a . 

Similarly the change in energy density and pressure 
through pair annihilation leads to a change in amplitude. 
In the WKB approximation, one finds that for kr a 1 
the amplitude C of $(fc, r) changes through pair annihi- 
lation by a factor 



C(k) 



a( T i) 
a(r 2 ) 





f 


'gesin)' 


1/2 




T2 


.5off(r 2 )_ 





= 0.911 



(62) 



where T\ <C t q <C t 2 . The WKB results for a(k) and 
C(k) match the numerical results shown in Fig. [7J 
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FIG. 8: Amplitude functions (solid curve) and f^ 3, 

(dashed curve) for the gravitational CDM density perturba- 
tion v of Eq. (I63L for r 2> T a , plotted versus wavenumber. 
For kr a < 1, / 3 pa = 1 and / 4 pa = C - \ = 0.0772 ■ ■ ■ . This 
figure assumes that the CDM is always collisionless and was 
never coupled to the radiation. The results show there is a 
7% enhancement of growth for waves that enter the horizon 
during pair annihilation. 



The change in the gravitational potential caused by 
pair annihilation induces changes in the CDM growth. 
The late-time solution of (|38ajl - l|38b|) with T/m x = is 



9C 



Hu 



cos 6 + 8 sin 9 



3C 



cosaCi(60 + sinaSi(6») 



+/s [hill — ) 

[smd-f 3 e] , 



(63) 



where 6 = 6 + a and the subscript indicates T/m x = 0. 
After pair annihilation the transfer functions fa and 
depend on k but not on r. 

Kinetic decoupling introduces a length scale r^. Pair 
annihilation introduces a second scale, r Q . When t q ^> 
these two physical effects can be separated. Thus we con- 
sider first the case of purely collisionless CDM, i.e. 7 = 
in Ij38(l , in which case there is no kinetic decoupling. Now 
there is only one timescale in the problem, r Q , and we de- 
note the corresponding transfer functions in i|63|) by /g a 
and / pa . The results obtained by numerical integration 
are shown in Fig. EI For kr a < 1, /f a « 1 - 3.00(fcr Q ) 4 
(where the coefficient 3.00 was found by numerical in- 
tegration) and fl a w C - i. For kr a > 1, the WKB 
approximation yields 



/•pa 



ft 



3 

pa 



1 - 0.1067(fcr a )- 



- 1 



ln(9.82fcr a ) , (64) 



where the coefficients 0.1067 and 9.82 were determined 
by numerical integration. The errors made by assuming 
are less than 1% for kr a > 100. 



For WIMP dark matter, the effects of kinetic decou- 
pling and pair annihilation combine to give transfer func- 
tions / 3 kd+pa and /4 d+pa in @3). For T/m x = and 
Td -C r a <C r these are 

/3 d+Pa (fc) = /i(fc)/ 3 Pa (fc) - 



kd+pa 



(k) 



C(k) 



h{k) -C + 



(65) 



[fi(k) 1] 



C(k) 



ln(1.589fer ) -lnC(fc) 



The coefficient 1.589 was found numerically but other- 
wise /3 d+pa and /4 d+pa were determined analytically. 
The general transfer function is given in terms of sepa- 
rate transfer functions for kinetic decoupling (/1, / 2 ) and 
pair annihilation (/ pa , f^)- 

Pair annihilation modifies the CDM velocity and den- 
sity transfer functions at k ~ r^ 1 due to the peak in 
/ pa (fc). There is also a small modification of the density 
for kr a ^> 1. In the limit fcr ^> kr a ^ 1, 
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h(k) 



In 61 + 0.08909 In 



7.416r a 



h{k) . 



(66) 

Comparing this with H40|) . we see that the term propor- 
tional to 1 — C(oo) = 0.08909 is induced by pair an- 
nihilation. This effect arises from the reduction of the 
gravitational potential by a factor C(oo) = a(ri) / cl(t2) 
at short wavelengths after pair annihilation. However, 
the correction to v is always less than 9%. 

The full density transfer function for nonzero CDM 
temperature during the radiation-dominated era follows 
from 1|48|) . (|63ll . and (|65ll . Numerical tables of the trans- 
fer functions /i(fcr d ), / 2 (fcr d ), C(fcr a ), a(kr a ), /£ a (fcT ), 
and / pa (fcT a ) are available from the author. 



V. EVOLUTION THROUGH 
RADIATION-MATTER EQUALITY 

Following kinetic decoupling and pair annihilation 
there is one additional major life event for CDM density 
fluctuations before they collapse to form nonlinear struc- 
tures: the transition to a matter-dominated universe oc- 
curring at l+z cq = a~q ss 3200. The constituents include 
photons, neutrinos, baryons, and CDM. The background 
equation of state is modified by the presence of nonrela- 
tivistic baryons and CDM, 



3w = 



1 



1 



y 



y = 



Pb + Pc 

(g cS /2)p~ l a, 



(67) 



eq 



Assuming that dark energy and spatial curvature can be 
neglected during the times of interest, the solution of the 
Fricdmann equation for t 3> r a is 



1 



r 

2r7 



o ff 2 



(68) 
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Here VL m is the present-day density parameter for p m = 
Pb + Pc- 

Before recombination, Thomson scattering couples the 
photons and baryons so that they behave as a single fluid 
with sound speed given by 



3c 



jb 



4p 7 



(69) 



We make two approximations which enable us to re- 
duce the dynamics to a fourth-order system in time, for 
which we obtain limiting analytic solutions: Photons and 
baryons behave like a single perfect fluid and neutrino 
shear stress is neglected. (We also assume a flat back- 
ground, K — 0, but this is always a good approximation 
during the times of interest.) These approximations in- 
troduce small errors in the results which may be elimi- 
nated by integrating the full system of equations for pho- 
tons, neutrinos, baryons, and CDM with CMBFAST 
or equivalent starting after pair annihilation from initial 
conditions obtained in the preceding section. However, 
the simplified dynamics leads to analytic results making 
it easy to distinguish the various physical effects in the 
transfer functions. 

With these assumptions and the further restrictions 
t r a and Td/m x <C 1, the perturbed Einstein equa- 
tions may be combined to yield an equation of motion 
for the gravitational potential $ = "J, 



$ + 3(l + cL)H$ + 3( 



2*6 - w)H 2 <S> + k 2 c 2 b $> 

PbPv j 



-A c c 2 b | S, 



PcP~i 



,(70) 



where S c and S v are, respectively, the CDM and neutrino 
density perturbation in conformal Newtonian gauge, and 
we have defined 



A c = inGa p c = - 



3 (l - h)y 



2 1 + y 



H 



(71) 



where fb = flb/Q m - Equation H70(l is a modified form of 
(JSBJ. The CDM fluid equations l|3"%|) remain valid with 
S = 5 C , u = u c . 

The evolution of metric perturbations is coupled to the 
evolution of both CDM and neutrino perturbations. To 
simplify the presentation we avoid solving the collision- 
less Boltzmann equation for neutrinos. Instead we make 
an additional approximation for the neutrino dynamics: 
either the neutrino fluid evolves like the photon-baryon 
fluid, or neutrino perturbations are damped by Hubble 
expansion on sub-horizon scales in a manner consistent 
with their evolution on super-horizon scales. In the first 
case we set S v — 5 7 on all scales, which is equivalent 
to setting S v = in Eq. I|70|l and replacing c 7 6 with c r b 
defined by 



3c rb — 



Pb 



4 (Scff/2)^ 



l + - 4 hy 



(72) 



In the second case, we note that on scales much larger 
than the Hubble length, for isentropic (curvature) per- 
turbations all species have the same number density per- 



turbation, S r 



5 V 



The evolution of isentropic 



perturbations on large scales follows from Ref. |14| : 
2 



S v 



1 



($ + ^- 1 $) . 



(73) 



On scales much smaller than the Hubble (or neutrino 
free-streaming) length, during the radiation-dominated 
era S v undergoes damped oscillations while the ampli- 
tude of the photon-baryon oscillations is constant. This 
qualitative behavior is correctly reproduced if one ap- 
plies (|73|l to the neutrinos at all times while leaving l|70|) 
unchanged. 

Thus, our two approximations for neutrino evolution 
correspond to two different choices for the sound-speed of 
the photon-baryon gas: either the sound speed is set to 
c rb with S v = or the sound speed is c 7 f, with S v given by 
(|73|l . Either way, when the CDM dynamics is added, the 
evolution reduces to a fourth-order system in time. The 
difference between the two treatments gives a measure of 
the importance of massless neutrinos for CDM evolution. 

Setting S v — and writing the sound speed as Cb which 
may be set to cither c 7 & or c rb , differentiating (|70|) and 
using the CDM continuity equation gives 

3, 



<9 t 3 $ + 5H$ + ^$[3 - 5io - 3(^(1 + w)] + 



3A c c 2 b <5> 



k 2 cl{<!> + H<£) = A c k 



4> "c 



(74) 



Differentiating this again and using the velocity equation 
u c + Ti-Uc — <& gives 



9 t 4 $ + (8 - 3c 2 b )Hd?<S> + 



17- 15(w + cg) 



(2-3cg)(l 



$ + 3W$ 



3») + T 



c rb 



I 2- 



wH 2 <S> 



Qwc 2 


n 2 <i> 


c 2 


c rb - 








4b)\ 


H 3 <f> 







= . (75) 



For wavelengths longer than the Hubble length this equa- 
tion is correct if c& = c r &; for wavelengths much shorter 
than the Hubble length it is approximately correct if 
Cb = c yb . We now present the analytic solution to the 
fourth-order system in these two limits. 

In the long- wavelength limit k 2 <§C TL 2 , the four linearly 
independent solutions may be obtained as quadratures 
using the methods of Ref. [l4| : 



4>4 



y 3 



y 2 (l + w) dr , 



y 

3H 

V 

n r 3 . 

-2 J y wdr , 

l~t f T 

-2 J y 2 w[l + {pb/p m )y} dr 



(76a) 
(76b) 
(76c) 
(76d) 
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These solutions give the time-dependence in the limit 
k — ► valid in both the radiation- and matter-dominated 
eras. The physical solution for isentropic curvature fluc- 
tuations is 



200 - 



100 



$4 



7*1 



9/2 8 16 
10 V + 9y"V~V 



1 



1 

16 



y for y <C 1 



(77) 



~5> 



-100 



-200 - 



Mode 1 is the decaying mode. Modes 3 and 4 are entropy 
perturbations and (for k 2 -C Tt 2 ) are constants added to 
the other solutions, 




2$, 



$4 



1 



" ftn 



(78) 



In the opposite limit, k 2 3> 7i 2 , approximate solu- 
tions to l|75|) may be found using the WKB method. The 
first-order WKB approximation gives two oscillatory so- 
lutions, 



FIG. 9: Results for the gravitational potential transfer func- 
tion at y = a/acq = 10 after the radiation-matter transition, 
assuming T d = 22.6 MeV, Q m h 2 = 0.131, Q b h 2 = 0.022. 
The solid black curve shows the effects of pair annihila- 



tion at T « 0.5 MeV as the bump at kr e 



10 5 - 8 and 



the effects of kinetic decoupling (with Td/m x = 0) as the 
damped oscillations at kr e > 10 7 ; the rapidly damped dotted 
curve includes the effects of free-streaming damping assuming 



, 2 3/2\_i f T / curve includes the effects of free-streaming damping assuming 

*± = (a c b ) expiifcy c b d,T , (79) (2T d /m x ) 1/2 = 0.0213. Other curves show the effect of elim- 



and two non-oscillatory solutions of the equation 
$ + 3W$ + I J- - 3 ) wH 2 $ = 

Z rb 



(80) 



This equation is equivalent to the usual evolution equa- 
tion for cold dark matter perturbations calculated assum- 
ing that the other components are unperturbed: 



inating either pair annihilation or kinetic decoupling. The 
straight dotted line is the WKB solution of Eq. I)86[l . verti- 
cally offset for clarity. 



gauge density perturbation v (not to be confused with 
the index of the Legendre functions) valid through the 
matter-radiation transition. On large scales, 



D + HD = A C D , D oc S cx y<$> 



(81) 



+ • 



2(fcr e )V$ 
4 + 3y 



For a(r) given by ijfilfy . the exact solutions of this equa- 
tion are 181 



where we have chosen the normalization $ = 
0. On scales shorter than the Jeans length, 



k 2 r 2 < 1 , (85) 
1 as fcr — > 



D + (t) =P v {l + r/2T e ) , TJ_(r) = Q„(l + r/2r e ) , 

V = - - + - \ ■!'■■> 24 1 1, . (82) 



i + iv/25-24/ b 



where -P„ and are Legendre functions of the first and 
second kind. Early in the radiation-dominated era they 
give 



9C/ 3 



D 



D, In 



Akr e 



V3 

cos(0 + a) 3(1 — /b)^ 



(kCbTe 



y 



2(kr e ) 2 y 



(86) 



1 - 



v{v + 1) r 



2 V r 



In the opposite limit, 



for t <C r e . 
(83) 



7rr(i/+ 1) 

1 \ T 



for r 3> T e 



(84) 



By matching the solutions found here to those in the 
radiation-dominated era, (|59|l and (|63|l . we obtain re- 
sults for the gravitational potential $ and synchronous 



where C = C(k) and /3 = /3 d+pa (fc) were given in the 
previous section. For k -C r^ 1 , C = / 3 = 1. 

Figure [5] shows the results obtained from numerical 
integration of (|?0"j) using lf73*j) to approximate the neu- 
trino evolution. In order to illustrate the logarithmic 
wavenumber-dependence of the growing mode in the 
matter-dominated era, the potential is scaled by — (fcr e ) 2 
and is plotted in a semilog fashion. If instead of using 
(|73|l the neutrinos are assumed to evolve like the photon- 
baryon fluid, the maximal change in $ is about 4% oc- 
curring at kr e rs 2 with negligible differences at much 
higher and lower frequencies. Thus the CDM transfer 
function is relatively insensitive to the detailed dynamics 
of massless neutrinos. For the CDM and baryon abun- 
dances used in the calculations, D + = 11.34 at y = 10. 
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The WKB result helps explain how gravitational po- 
tential fluctuations with an rms amplitude 2 x 10~ 5 on 
large scales can generate nonlinear structures at high 
redshift. Growth after the universe becomes matter- 
dominated contributes a factor D + 1 + %y once the 
baryons are released and join the CDM potential wells. 
Evolution after kinetic decoupling in the early universe 
contributes an additional factor 91n(4feT e /-\/3)- The 
logarithm is sometimes called the Meszaros effect [Tij : 
it arises because acoustic oscillations in the relativistic 
plasma gravitationally induce velocity perturbations in 
the dark matter. Half of the factor of 9 (i.e., a factor of 3) 
arises because the comoving number density is affected 
by the cube of the strain (1 — ty) in three dimensions. 
Altogether, CDM density perturbations in the matter- 
dominated era are enhanced by a factor ^yln(4fer e /\/3) 
which is sufficient to drive rms density perturbations non- 
linear on scales below 10~ 5 M Q by a redshift of 20. 

These results neglect free-streaming of the CDM par- 
ticles. As shown in Section IIIII during the radiation- 
dominated era CDM free-streaming leads to a modified 
Gaussian suppression, Eq. (|48|) . Including the radiation- 
matter transition, this gives 

(tOT) - /|Wl±fi^), (87) 

W 5m x V 1 + 4r e /T J 

where r* w r^. For the parameters of Fig. [5] with pair 
annihilation, r e /rd = 3.895 x 10 7 . 

VI. SMALLEST-SCALE CDM STRUCTURES 

Having computed the evolution of the CDM transfer 
function through kinetic decoupling, pair annihilation, 
and radiation-matter equality, we can now make pre- 
dictions for the smallest-scale CDM structures. To be 
definite, we adopt the standard flat ACDM model with 
tt b h 2 = 0.022, Q c h 2 = 0.109, h = 0.71, f2 A = 0.740. 
With these assumptions, a^ 1 = 3160, r e = 147 Mpc, 
Td = 3.78 pc, and k^ 1 = 1.18 pc. 

The CDM density transfer function v(k,r) is normal- 
ized so that the variance of 5p/p in the usual synchronous 
gauge is a 2 = J v 2 {k, r)Pi(k) d 3 k where Pi(k) is the spec- 
tral density of the gravitational potential $ early in the 
radiation-dominated era. The radiation-era potential is 
related to the curvature perturbation in a flat universe 
C = K = k by $ = §«. 

The treatment of this paper makes several approxima- 
tions: neutrino shear stress is neglected; the Boltzmann 
equations for neutrinos and photons after recombination 
have not been integrated; and photons and baryons are 
assumed to be tightly coupled until y = a/a cq = 10 
(z = 315). To test this crude treatment of multiple 
fluids, we compared the resulting transfer function with 
that computed numerically from CMBFAST 0. For 
kr e w 100 (small enough so that pair annihilation and 
kinetic decoupling are unimportant), our v(k) was too 
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1 ' 1 1 ' 1 ' L - 

-8 -6 -4 -2 2 4 

log 10 M/M sun 



FIG. 10: RMS mass density perturbation in a sphere con- 
taining mass M, at redshift z = 30.6 (a/a cq — 100). Two 
different assumptions are made for the scalar spectral slope 
n 3 . For each case, three different choices are shown for 
the neutralino thermal speed at decoupling: (2Td/m x ) 1//2 £ 
{0,0.01,0.0213}. 

small by 7% (due to the delay in baryons joining CDM 
until z = 315) while for fcr e < 0.1 our v(k) was too large 
by 3.6% (due to the neglect of neutrino and photon shear 
stress). A more accurate calculation would incorporate 
our transfer functions for pair annihilation and kinetic 
decoupling into CMBFAST or equivalent code, but this 
level of accuracy is unnecessary for the present purposes. 

Figure ITU1 shows the rms mass density perturbation in 
a sphere containing mean mass M, which follows from 

a 2 {M) = J v 2 {k,T)W 2 {kR)P l {k)d i k , (88) 

where W(x) = 3j±(x)/x is the spherical tophat win- 
dow function and R = (3M/4np m ) 1/3 - The primordial 
spectrum has been normalized to A 2 ^ = 2.40 x 10 -9 at 
k = 0.002/Mpc from the WMAP 3-year results [lr| . 
The amplitude for the scale-invariant case n s = 1 is 
more than 50% greater than it is in the tilted model 
n s = 0.94 favored by the data. With a neutralino mass of 
100 GeV and kinetic decoupling temperature 22.6 MeV, 
(2T d /m x ) 1 / 2 = 0.0213. In this case, 3(7 perturbations 
reach v = 1 by z = 47. 

The Press-Schechter mass function dn/dM obeys 
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(2\ 1/2 
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d\nM 
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dhxM 





where 8 C = 1.686. The results, shown in Fig. 1111 suggest 
that the smallest CDM objects may have masses even 
smaller than an earth mass, 3.0 x 10~ 6 Mq. 

WIMP decoupling imprints two different length scales 
on the spectrum: the comoving horizon size oc T^ 1 
and the comoving free-streaming distance (?d/m x ) 1 / 2 T ci ; 
the former scale is larger and more important. The char- 
acteristic mass for the smallest objects is the CDM mass 
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log 10 M/M sun log 10 M/M, 



FIG. 11: Press-Schechter mass fraction for the cases shown in 
Fig. 1101 The effects of increasing free-streaming damping are 
apparent in the erasure of the smallest objects with increasing 
(2Td/m x ) 1 ^ 2 . The cases of no free-streaming damping have 
bimodal distributions, however this requires unreasonable as- 
sumptions about the dark matter mass and couplings. 

in the horizon at decoupling, 
M d = —p m {a = l)(cT rf ) 3 

= 7 . 59 * 10 -^(^)" 5,4 M S , < 9 0) 

where C was defined in With the default values 

C = 0.0433 and g cB = 10.75, M d = 8.29 x 10~ 6 M , or 

M A = 1.0 x 10~ 4 ( — — — 1 M ra , (91) 
V 10 MeV / ° ' v ' 

in agreement with Ref. . Berezinsky et al. [f| derived a 
different characteristic mass assuming that diffusion fol- 
lowed by free-streaming damping sets the minimum ob- 
ject mass. As we have shown, for {2Td/m x ) 1 ^ 2 <§; 1, the 
dominant damping process is not diffusion or free stream- 
ing but friction between WIMPs and relativistic leptons 
during kinetic decoupling (Silk damping). That this is 
not a diffusive process may be seen by its persistence in 
the limit Td/m x = 0. 

The evolution of the mass function requires comput- 
ing the density perturbations at smaller redshift. The 
baryons present a slight complication because after re- 
combination, even though they are free to move apart 
from the photon gas, they resist gravitational instabili- 
ties for k > fcjb where the baryon Jeans wavenumber is 

fc.j b =(^y /2 ^, 02) 

V 2a J Cb 



FIG. 12: Press-Schechter mass fraction for n s — 0.94, m x = 
100 GeV and T d = 22.6 MeV, for several redshifts. Con- 
stant (M / p m )dn/d In M indicates equal mass is contained in 
objects whose masses span equal logarithmic mass intervals. 
The slight bump at 50 Mq is due to enhanced growth during 
electron pair annihilation. The growth on scales M > 10 J 
Mq has been underestimated by ignoring the boost given by 
infalling baryons. 



with Cb being the baryon sound speed. The small residual 
ionization fraction persisting after the nominal recombi- 
nation maintains the baryon temperature close to the 
microwave background temperature until a « axdec = 
1/125 yielding (for the standard cosmological param- 
eters) kjh — 252 Mpc -1 for a <C aTdcc and /cjb = 
252-\/a/aTdec Mpc -1 for a ^> aTdcc • Thus, even before 
rcionization, baryon perturbation growth is inhibited on 
scales k > 252 Mpc -1 or dark matter masses less than 
about 10 5 M-q. Since we are mainly interested here in 
smaller scales, we will assume that baryons do not par- 
ticipate in gravitational instability. The CDM density 
perturbations therefore grow with time according to the 
D+(t) solution JgSJ with f b = Q b /Q m . For f b = 0.168, 
at a = 1 this results in an almost 50% reduction in lin- 
ear growth compared with f b = 0. Dark energy plays a 
role only at very late times. At z = 2.16, a cosmological 
constant with J7a = 0.74 suppresses the growth of linear 
density perturbations by 2%. 

Figure^|shows the evolution of the mass fraction with 
redshift for plausible WIMP parameters. At z = 31 
the peak of Mdn/dlnM occurs at M w 2.3Md (6 
earth masses) but there may be numerous smaller ob- 
jects. Indeed, the numerical results imply da/d\nM oc 
-(M/M d ) 2/3 for M < M d . This follows from W(kR) 
being an analytic function of (kR) 2 with W(0) = 1. For 
R = the variance integral (|88|l converges because of the 
small-scale damping caused by kinetic decoupling. Thus, 
analyticity implies a 2 oc C\ — CiiRjTd) 2 as R — > where 
C\ and C2 are independent of filter radius R oc M 1 / 3 , 
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hence da/dlnM oc —(R/rd) 2 as R — > 0, or 



dn 



d\nM 



oc M 



-1/3 



(93) 



Given a smooth, symmetric window function, the number 
of small objects diverges but the mass contained in them 
converges. 

These results for M < Md are uncertain because there 
are window functions W(kR) that give a different result. 
For example, a sharp fc-space filter, with W(kR) = 1 for 
kR < 7r and zero otherwise, when combined with the free- 
streaming damping of v(k, r), leads to an exponential 
cutoff instead of a power-law rise of dn/d In M as M — > 0. 
However, a sharp fc-space window gives rise to long-range 
oscillations in physical space, which seems unlikely to 
properly describe the local physics of dark matter halo 
formation. 

The choice of window function represents an enormous 
simplification of the nonlinear dynamics of dark matter 
halo formation. N-body simulations are needed to deter- 
mine the correct mass function for M < Mj,. Such simu- 
lations must fully resolve scales below the free-streaming 
damping length of 1.17 pc in the initial mass distribu- 
tion. Cold dark matter caustics 0] are also likely to be 
important in the formation and dynamics of the smallest 
halos. 

By redshift 9, RMS density perturbations become non- 
linear and the small-scale mass distribution is nearly 
scale-invariant over 8 orders of magnitude of mass, 1CP 5 
M Q to 10 3 M & , reflecting the nearly scale-invariant (n w 
—3) character of the small-scale CDM density field. At 
later times larger objects build up by mergers of sub- 
structures as the effective spectral slope n > — 3 on larger 
scales. The smallest-scale structure forms by almost sud- 
den collapse as opposed to hierarchical clustering. 

The results presented here do not include tidal destruc- 
tion of small objects formed at high redshift. The sur- 
vivability of earth-mass halos remains an open question 
0, El> E3 ■ The vast dynamic range of masses reflected 
in Fig. 1121 provides an enormous challenge to numerical 
simulation methods atte mpti ng to determine the survival 
of the first forming halos [21j ■ 



a solar mass, with precision similar to cosmic microwave 
background measurements on degree angular scales, we 
could hope to determine the WIMP mass and elastic scat- 
tering cross section with leptons by astrophysical means. 
This seems far-fetched because of the tremendous degree 
of nonlinear processing that has occurred since the first 
cold dark matter structures formed at a redshift about 
50. However, if any of the smallest bound objects survive 
tidal stripping, theymay produce an observable WIMP 
annihilation signal |7| which provides a window into par- 
ticle physics. 

Moreover, the nearly scale-invariant structure pre- 
dicted for dark matter fluctuations on mass scales below 
10 4 M Q could affect the formation of larger-mass struc- 
tures. This mass range is very difficult to study with 
standard numerical simulation methods. 

While the present work has answered the question of 
how small-scale fluctuations in the WIMP distribution 
evolve during the linear stage of evolution, it raises new 
questions about nonlinear halo formation. If the Press- 
Schechter formalism with a local spatial filter (or analytic 
window function) correctly describes the formation of the 
smallest halos, then there is no minimum mass for WIMP 
microhalos; the mass function dn/dlnM oc Af" 1 / 3 for 
M < Kr 4 (T d /10 MeV)~ 3 M Q . Most of the mass, how- 
ever, is in larger clumps. Nonlinear evolution is likely 
to strongly modify the Press-Schechter distribution. Not 
only is tidal stripping important, but the Press-Schechter 
model is based on the assumption of hierarchical cluster- 
ing, which breaks down on the small scales considered 
here. 

The present paper has considered dark matter only 
in the form of a thermal relic, i.e. particles that were 
once in thermal equilibrium with the relativistic plasma 
of the early universe. The dark matter may be instead 
a nonthermal relic — the axion — which began its life 
as a Bose-Einstein condensate. The axion field begins 
oscillating after the QCD phase transition but is never 
collisionally coupled to the plasma. It would be interest- 
ing to investigate the smallest-scale structure in an axion 
dark matter model. 
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